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A minimal model for excitons within time-dependent density-functional 
theory 
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The accurate description of the optical spectra of insulators and semiconductors remains an important chal- 
lenge for time-dependent density- functional theory (TDDFT) . Evidence has been given in the literature that 
TDDFT can produce bound as well as continuum excitons for specific systems, but there are still many 
unresolved basic questions concerning the role of dynamical exchange and correlation (xc). In particular, the 
role of the long spatial range and the frequency dependence of the xc kernel / xc for excitonic binding are 
still not very well explored. We present a minimal model for excitons in TDDFT, consisting of two bands 
from a one-dimensional Kronig-Pcnncy model and simple approximate xc kernels, which allows us to address 
these questions in a transparent manner. Depending on the system, it is found that adiabatic xc kernels can 
produce a single bound exciton, and sometimes two bound excitons, where the long spatial range of f xc is not 
a necessary condition. It is shown how the Wannier model, featuring an effective electron-hole interaction, 
emerges from TDDFT. The collective, many-body nature of excitons is explicitly demonstrated. 

PACS numbers: 31.15.ee, 73.20.Mf 



I. INTRODUCTION 

The study of the electronic structure of materials usu- 
ally begins with noninteracting electrons due to the vast 
number of particles involved. Many-body methods then 
provide a hierarchy of corrections to account for the 
Coulomb interaction to various order. Many-body ap- 
proaches such as GW 1 ' 2 and the Bethe-Salpeter equation 
(BSE) 3,4 are frequently and successfully employed in the 
calculation of the electronic structure and excitations of 
materials. Though accurate and physically sound, these 
many-body methods can become cumbersome and im- 
practical for large systems due to the steep scaling of the 
numerical cost versus the system size. 

Alternatively, density-functional theory (DFT) and 
time-dependent density- functional theory (TDDFT) 
are popular methods for calculating electronic ground 
states and excitations, respectively, and are widely used 
in chemistry, physics, materials science, and other ar- 
eas. Density-functional methods solve the many-body 
problem by constructing a noninteracting system which 
reproduces the electronic density of the interacting, phys- 
ical system. The favorable balance between accuracy and 
efficiency makes the resulting DFT and TDDFT schemes 
unrivaled for large but finite system sizes. 8 

Considerable effort has been spent to replicate this 
success of TDDFT for periodic solids. 9 Generally speak- 
ing, TDDFT works very well for simple metallic systems, 
where the excitation spectrum is dominated by collective 
plasmon modes. The reason is that common local and 
semilocal exchange-correlation (xc) functionals are based 
on the homogeneous electron liquid as reference system, 
which is an ideal starting point to describe electrons in 
metals. 

The situation is more complicated in insulators and 
semiconductors. The first problem that comes to mind 
is that of the band gap, which is typically strongly un- 



derestimated by most popular xc functionals of DFT. In 
principle, TDDFT provides a mechanism to obtain the 
correct band gap, 10-12 but this puts very strong demands 
on the xc kernel / xc (it has to simulate a discontinuity, 
which requires a strong frequency dependence). 

The second difficulty are excitonic effects. It is a well- 
known fact that standard local and semilocal xc func- 
tionals do not produce any excitonic binding; 4,9 again, 
the proper choice of / xc is crucial. There are many exam- 
ples in the literature of successful TDDFT calculations 
of excitonic effects, using exact exchange, 13 an effective 
xc kernel engineered from the BSE, 14-18 a meta-GGA 
kernel, 19 and a recent 'bootstrap' xc kernel. 20 These ker- 
nels all have in common that they have a long spatial 
range; however, it has also been shown that certain exci- 
tonic features can be equally well reproduced by simple 
short-range kernels. 9,15 This calls for further explanation. 

Due to the complexity of real solids, the question 
of the general requirements for excitonic binding in 
TDDFT has been difficult to analyze. As a first step 
towards a simplified TDDFT approach for excitons, a 
two-band model was recently developed, which was used 
to test the performance of simple xc kernels for calcu- 
lating excitonic binding energies in several III-V and II- 
VI semiconductors. 21 ' 22 In this paper we will push this 
reductionist approach further and propose a minimal 
TDDFT model for excitons. 

Our model is one-dimensional (ID) and uses two sim- 
ple Kronig-Penney-type bands as input. We show that 
the minimal model reproduces and reveals many aspects 
of excitons. The model is accessible and relatively easy 
to implement, and it can be used to identify impor- 
tant aspects of the xc functional for excitonic effects. It 
clearly shows that excitons are collective excitations of 
the many-body system: the appropriate phase-coherent 
mixing of the single-particle excitations is accomplished 
via a coupling matrix featuring / xc . The properties of 
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FIG. 1. Schematic optical spectrum of a typical 3D direct-gap 
insulator. Dashed line: independent-particle spectrum. Solid 
lines: spectrum including excitonic effects. 



this coupling matrix are analyzed and compared with its 
BSE counterpart. 

In textbooks, excitons are usually introduced through 
the two-body Wannier equation, which describes an elec- 
tron and a hole which interact via a screened Coulomb 
potential. While this arguably constitutes the simplest 
model for excitons, it is based on several drastic assump- 
tions which are not fulfilled in general. 23 We discuss how 
and under what circumstances a Wannicr-likc equation 
emerges from our minimal TDDFT model. 

The paper is structured as follows. We give an intro- 
duction to Wannier excitons in Sec. II A and to TDDFT 
for solids in Sec. II B, followed by a description of the 
minimal model in Sec. II C. Section III then presents re- 
sults for the minimal exciton model comparing TDDFT 
with the BSE, and discusses various implications. In 
Sec. IV we show how the Wannier equation emerges from 
TDDFT. Conclusions are given in in Sec. V. Details of 
the BSE method are provided in Appendix A. Atomic 
units {h = e = m e = (47reo) 1 = 1) are used through- 
out unless otherwise stated, and we will only consider 
spin-unpolarized systems. 



II. BACKGROUND AND MODEL 
A. Wannier excitons 

The electronic structure of crystalline solids is de- 
scribed by the Bloch theory, where the electrons move in 
a periodic effective single-particle potential which reflects 
the crystal symmetry. As a result, the electronic states 
form energy bands. In insulators and semiconductors, 
electronic excitations take place between the occupied 
(valence) and unoccupied (conduction) bands. These 
interband transitions can be described within a simple 
independent-particle approach based on Fermi's Golden 
Rule; one thus obtains a reasonable qualitative account 
of the optical properties in these materials. 24 

However, experiments reveal that there are important 
modifications to this picture, as illustrated schematically 
in Fig. 1. Above the band gap, the spectrum appears 
strongly enhanced, and below the band gap one may find 




FIG. 2. Illustration of Wannier exciton as an electron-hole 
pair, extending over many lattice constants. 

discrete absorption peaks known as bound excitons. The 
origin of these modifications are Coulomb interactions: 
the simple picture of independent single-particle excita- 
tions is replaced by a more complex scenario where these 
excitations are dynamically coupled. Excitonic effects 
are ubiquitous in nature, and occur in 3D, 2D and ID 
systems alike. 25 ' 26 The details of the excitonic modifi- 
cations to the noninteracting spectrum depend on the 
dimensionality of the system. 27 

The standard textbook explanation of excitons is based 
on the simple picture of an electron-hole pair held to- 
gether by Coulomb interactions, see Fig. 2. One thus ar- 
rives at a two-body problem similar to the positronium 
atom, with a center-of-mass momentum k and relative 
motion described by the Wannier equation: 27 



2m, 



V(t) 



Vv(r) = E u ip v (r), 



(1) 



where m r is the reduced mass, defined as m" 1 = m~ — 
m~ l . m c and m v are the effective masses of the conduc- 
tion band electrons and the valence band holes. V(r) = 
1 jer is the Coulomb interaction between the electron and 
the hole, divided by the static dielectric constant of the 
system. i\) v and E v are the excitonic wave function and 
binding energy, respectively. 

In ID systems, the Coulomb interaction is ill-defined 
and requires, in general, some parametrized form; 28 we 
will use the following soft-Coulomb interaction: 



A 



(2) 



where A and a are parameters. In the following, we will 
use a = 0.01 throughout. 

The Wannier equation (1) has the form of a hydro- 
genic Schrodingcr equation, so it possesses a Rydberg 
series (even for ID cases with soft- Coulomb interaction) 
with infinitely many eigenvalues below the band gap, and 
a continuum above the gap. 27 The excitons in the Wan- 
nier model for 3D and 2D cases both enhance the optical 
absorption near the band gap, while the presence of exci- 
tons in ID systems suppresses the optical absorption just 
above the band gap. 



3 



Sham and Rice 23 showed how the Wannier equation 
can be derived from first principles starting from the 
BSE, under the assumption that the Bohr radius of the 
exciton is much bigger than the lattice constant. The 
resulting Wannier picture of excitons appears clear and 
intuitive, but this simplicity is somewhat deceptive. In 
reality, excitons are a dynamical many-body phenomenon 
and require a subtle coordination and cooperation of 
many single-particle transitions between two bands. In 
the following, we will develop a model based on TDDFT 
which will illustrate the true physical nature of excitons, 
but which will remain sufficiently transparent to allow a 
simple interpretation of the collective many-body effects 
that are responsible for the excitonic binding. 



energies, i<Hxc m Eq. 



? (ij)(mn) 



is defined as 



d 3 r / dV &(r)^(r)/ HX c(r,r» 

x C(r')<Mr'), (6) 



where the 4>'s are the ground-state Kohn-Sham orbitals, 
and /hxc is the Hartrce-exchange-corrclation (Hxc) ker- 
nel, defined as a Fourier transform of 
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B. TDDFT in finite and periodic systems 

TDDFT 7 is an in principle exact approach for elec- 
tron dynamics, based on the uniqueness of the mapping 
between the time-dependent electronic density and the 
external potential. 5 The key equation of TDDFT is the 
time-dependent Kohn-Sham (TDKS) equation: 



(3) 



where <pi are the TDKS orbitals of the nonintcracting 
Kohn-Sham system which reproduces the density of the 
real interacting system. w ext and i>h are the external po- 
tential of the physical system and the Hartree potential, 
respectively, and the xc potential v xc is the only piece 
that needs to be approximated in practice. 

The excitation spectrum of a system can be calculated 
via time propagation of Eq. (3) following a suitably cho- 
sen initial perturbation. Alternatively, one can obtain ex- 
citation energies and optical spectra directly from linear- 
response TDDFT, using the so-called Casida equation: 29 



The xc kernel / xc has to be approximated in practice. 

X and Y make up the eigenvector in Eq. (4), and 
they describe excitations and de-excitations, respectively. 
A commonly used approximation to Eq. (4), known as 
Tamm-Dancoff approximation (TDA), 30 is to set B = 0, 
so that the Casida equation reduces to 

\&i m 5j n (e J - £i) + F^ mn \uj) p mn (w) = u>pij(u). 

(jnn) 

(8) 

This decouples excitations and de-excitations, and the 
computational cost is reduced. There are situations, 
for instance for molecular excitations of open-shell sys- 
tems, 31 in which the TDA is preferred in practice over 
the full Casida equation (4). We find that the TDA can 
also be advantageous for excitons (see Sec. III). 

By considering only a single Kohn-Sham transition in 
Eq. (4) one arrives at the small-matrix approximation 
(SMA): 32 < 33 



J SMA,ij 



J KSaj 



(9) 



where loks ■ 



e, is the Kohn-Sham transition fre- 



quency to be corrected. One can further simplify this by 
making the TDA, which yields the single-pole approxi- 
mation (SPA): 



A B 
B A 



1 
1 



(4) 



Equation (4) is a generalized eigenvalue equation. One 
obtains the optical transition frequencies from the eigen- 
values (J, and the corresponding eigenvectors tell us how 
the Kohn-Sham single-particle transitions are mixed to 
form the transitions of the interacting system. The ma- 
trices A and B in Eq. (4) are defined as 



A MW (w) = Wjn(£ 3 -ei)+^ 
Bfc')(-")( tl ;) = F HXC 



? (ij){mn) 



(ij)(mn 

HXC 



(5) 



where i,j,m,n arc labels for ground-state Kohn-Sham or- 
bitals, and the e's are the associated Kohn-Sham orbital 



WSPA,i 



(10) 



The SMA and SPA are valid when the considered exci- 
tation is far away from other transitions in the system. 
Though they are not usually accurate enough for real 
calculations, their simplicity makes them very useful for 
theoretical analysis and development. 

In periodic solids, the Kohn-Sham orbitals are labeled 
with the band index i and the wavevector k and have the 
Bloch form: 



(r) 



e ikr M lk (r). 



(11) 



It is in principle possible to adapt the Casida equation (4) 
for the case of orbitals of the form (ll), 34 ' 35 and use this 
to calculate the excitation spectrum. However, to obtain 
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the optical spectrum in solids it is more convenient to 
calculate the macroscopic dielectric function: 4 



emM = lim 



1 



9^0 e 



G=G' 



= (q> w ) 
1 



(12) 



lim 

<?^o i + i) G=0 (q)xG=G'=o(q, w) 



where the G's are reciprocal lattice vectors, and x = 
dn/dv ex t is the linear response function, ejvf(w) can be 
expressed as 4 



cm(w) = 1 - lim «G=o(q) ^ 

q— >0 z — 4 



|i> 



uj - uj\ + ir) 



(13) 



where A labels the solutions of Eq. (8) for an extended 
system. 

Beyond linear response, few real-time TDDFT calcu- 
lations exist for periodic solids. 36,37 Instead of directly 
solving the TDKS equation (3), the TDKS orbitals can 
be expanded in terms of ground-state Kohn-Sham Bloch 
functions as 



ik(r,t) = y^Qmkft)0mk(r). 



(14) 



The time-dependent density matrix is defined as 

PTW = c mk (t)c* nk (t). (15) 
The equation of motion of the density matrix is then 

i^Pik(t) = [H k (t),p jU (t)], (16) 
with the TDKS Hamiltonian matrix Hk(£) defined by 

Hr(t)= f d 3 rC lk (r) J ff K s(r,O0nk(r), (17) 

where H is the volume of the unit cell, and HKs( r it) is 
the TDKS Hamiltonian of Eq. (3). This density-matrix 
approach has been used to derive the TDDFT version of 
the semiconductor Bloch equations. 21,22 In this formal- 
ism we only consider vertical transitions, where the Bloch 
wavevector k does not change during the dynamics. Non- 
vertical excitations are not considered, since they involve 
indirect (i.e., phonon-assisted) processes which we ignore 
here. 



C. Minimal model for excitons 

Solids are formally described by the many-body Schro- 
dinger equation. Since exact solutions are not possible, 
it is instructive to resort to model systems to eliminate 
undesired details of the many-body system and provide 
clear illustrations of the specific features one is interested 
in. The Wannier equation for excitons presented in Sec. 
II A is such a model. 



While intuitive, the Wannier model assumes the elec- 
tron-hole interaction as given, so the excitonic effects are 
already built in by default. However, this does not ex- 
plain under what conditions one expects to see the forma- 
tion of excitons, and the many-body nature of excitonic 
effects remains hidden. We therefore propose a minimal 
model for excitons which lowers the abstraction level of 
the Wannier model, and where excitonic effects show up 
without any ad hoc assumptions. 

For excitations near the band gap, a reasonable ap- 
proximation is to use a two-band model, i.e., only to 
consider the highest valence band (v) and the lowest con- 
duction band (c). This means that we only need to con- 
sider those elements of the time-dependent density ma- 
trix p^™(£) [Eq. (15)] for which mn = cc,cv,vc,vv and 
i = v; the latter index will be dropped in the following. 

For the case when a small perturbative electric field is 
applied to the system, it is sufficient, to lowest order in 
the perturbation, to consider only the time evolution of 
the off-diagonal part of the density matrix, 22 One 
then obtains from Eq. (16) 



<^HXC,k(*), 



,cv 

J k 



e c k - euk, and SV££ ck (t) 



yev 

r HXC.k 



(18) 

(*)- 



yev 

HXC.k 



(0). Here, V££ c k denotes the matrix elements of 



vn(r,t) +v xc (r,t), defined similarly to Eq. (17). There 
is no external perturbation in Eq. (18); we assume free 
propagation of the system. Our reasoning is that indi- 
vidual excitations can be viewed as the eigenmodes of 
the system and do not depend on the specific form of the 
perturbative field. Fourier transformation of Eq. (18) 
gives 



1 Z^k' \^HXC,k,k'Pk' l W J + ^HXC,k,k'Pk' ) 



where the factor 2 accounts for the spin, and 



7 (ij)(mn) 
HXC.k. k' 



d 3 r / dV ik (r)0* k (r) 
n Jn 

x / HXC (r,r',a;)0; ik ,(r')0 n k'(i-') 

((ijlfuxalmn)). 



(19) 



(20) 



Equation (19) is equivalent to the SMA for finite sys- 
tems. While the SMA only refers to the transition be- 
tween one individual occupied and one individual unoc- 
cupied orbital, Eq. (19) considers the transitions between 
the valence and the conduction band as a whole. Ignor- 
ing the coupling between excitations and de-excitations 
by setting F^"^ cv ^ = (i.e., making the TDA), one ar- 
rives at the solid-state analog of the SPA: 



E 

k' 



^Sm1^h=^h. (2i) 



This is the central equation which we will use to describe 
excitonic effects. It requires as input the Kohn-Sham 
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FIG. 3. ID Kronig-Penney model: 
and band structure (upper panel). 



potential (lower panel) 



Bloch functions for the valence and the conduction band 
of an insulator or semiconductor. To keep things as sim- 
ple as possible, we will consider the Kronig-Penney (KP) 
model 38 rather than a real material. The KP model is 
a ID nonintcracting system with a periodic potential of 
square wells. Within the unit cell [—6, a], the potential is 



V KP (x) 



< x < a 
V -b < x < 



(22) 



A typical example for the band structure of the KP 
model is plotted in Fig. 3. Despite its simple appear- 
ance, the KP model is very versatile: by varying the 
values of the lattice constant a + b, barrier width b, and 
barrier height Vq, a wide range of band gaps and band 
curvatures can be achieved. A square- well potential of fi- 
nite depth does not support an infinite number of bound 
states as the Coulomb potential does; but this, in fact, 
closely reflects the reality of the effective potential felt by 
the valence electrons in a solid, which is relatively shallow 
due to the screening of the bare nuclear charges by the 
core electrons. In practice, many solid-state calculations 
account for this screening by using muffin-tin potentials 
or pseudopotentials. 39-41 The KP model can be viewed 
as an elementary version of this approach. 

In the following, we always choose the first two bands 
to be fully occupied and the higher bands to be empty. 
We then make the two-band approximation for band 2 
and 3, since the shapes of these bands resemble the high- 
est valence band and lowest conduction band in direct- 
gap materials such as GaAs. The bands in the KP model 
are sufficiently well separated, so the two-band approxi- 
mation is justified. 

To establish a connection with TDDFT, we assume the 
solution to the noninteracting KP model as our ground- 



state Kohn-Sham system. In other words, the poten- 
tial in Eq. (22) represents the exact Kohn-Sham poten- 
tial v ex t + wh + v xc which corresponds to a physical sys- 
tem whose external potential w cxt is uniquely determined 
thanks to the Hohenberg-Kohn theorem of DFT. For our 
purpose, it is not necessary to know what this external 
potential looks like. The Kohn-Sham Bloch functions can 
then be determined in an elementary fashion. 38 

For comparison, we will also carry out BSE calcula- 
tions in our model system (see Appendix A for tech- 
nical details). BSE calculations are typically based on 
ground-state quasiparticle states obtained from the GW 
method. 4 This is because the single-particle gap in GW is 
usually closer to experiment than the approximate Kohn- 
Sham gap. However, in our case this distinction is not 
important because we use the given KP band structure 
as input for both BSE and TDDFT. 

In summary, our minimal TDDFT model for excitons 
consists of the following two ingredients: 

(1) A two-band model for the vertical transitions be- 
tween the highest valence band and the lowest conduction 
band, see Eq. (21); 

(2) the band structure from a ID KP model. 

Of course, the model is not complete without a choice 
for the xc kernel / xc . This will be discussed below. 



III. RESULTS FROM THE MINIMAL MODEL 

A. Bound excitons from the BSE and from TDDFT 

The exact xc kernel / X c( r j r'lW) is unknown and must 
be approximated; we restrict ourselves to adiabatic ker- 
nels that have no frequency dependence. The adia- 
batic local-density approximation (ALDA), as well as all 
scmilocal, gradient-corrected xc kernels, are known to be 
unable to describe excitonic effects. 4 The exact xc kernel 
has a long-range decay of 1/ |r — r'|, which is absent in 
all (semi)local xc kernels derived from the uniform elec- 
tron gas. This long-range part is thought to be essential 
for excitons. 4,7 ' 17 

The long-range behavior of f xc depends on the dimen- 
sionality, and in our ID model system we define the fol- 
lowing long-ranged (or 'soft- Coulomb') xc kernel: 



A sc 



y/(x ~ X 1 ) 2 + , 



(23) 



We also consider an extremely short-ranged contact xc 
kernel: 



f™ nt (x,x',uj) = -A cont S(x - x 1 ). 



(24) 



These model xc kernels depend on the constants A sc and 
y4. cont , which we will treat as fitting parameters in the fol- 
lowing. The idea is to tune the parameters in the model 
xc kernels so that bound excitons are produced, and to 
align the lowest bound exciton in the TDDFT spectrum 
with the lowest bound exciton in the BSE spectrum. 
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FIG. 4. Imaginary part of the dielectric function, Im(ejv/), 
calculated with BSE and TDDFT. Parameters of the KP 
model: a = 2.6, b = 0.4, V = 8. For BSE, A = 0.25 in 
Eq. (2). For TDDFT with the contact kernel, Eq. (24), 
A cont = 2.32. For TDDFT with the soft-Coulomb kernel, 
Eq. (23), A sc = 0.898. The BSE produces several bound 
excitons, but TDDFT only one. 



FIG. 5. Same as Fig. 4, for KP model parameters a = 3, b — 
3, V = 1. For BSE, A = 0.14 in Eq. (2). For TDDFT with 
the contact kernel, Eq. (24), A cont = 3.77. For TDDFT with 
the soft-Coulomb kernel, Eq. (23), A sc = 0.955. TDDFT 
produces two bound excitons. Higher-lying bound excitons 
exist within BSE but are numerically hard to resolve. 



Results for the imaginary part of the dielectric func- 
tion are presented in Figs. 4 and 5. We find that both 
the long-ranged an d the short-ranged f^° nt produce 
bound excitons, and thus, strictly speaking, the long- 
range behavior of the xc kernel is not really required for 
excitonic effects. The BSE results in Figs. 4 and 5 show 
several identifiable bound excitons, in agreement with 
the Rydberg series predicted by the Wannier model; the 
number of visible bound excitons somewhat depends on 
the numerical resolution in momentum space. 

For the KP model parameters of Fig. 4, we find that 
the adiabatic TDDFT can only bind a single excitonic 
state. For other KP parameters (specifically those in 
which the lowest conduction band is above the barrier), 
TDDFT produces two excitons, see Fig. 5, which agree 
well with the lowest two excitons in BSE. There are ad- 
ditional, higher-lying bound excitons in BSE which are 
very faint and difficult to resolve numerically. For all the 
KP systems we tested, we never found more than two 
bound excitons with TDDFT. This indicates the limita- 
tions of the adiabatic xc kernels used here. 

As mentioned in Sec. II C, the Wannier model does 
not clearly demonstrate that excitons are collective ex- 
citations. Since the Wannier model assumes a single 
electron-hole pair picture, one cannot immediately see 
that excitonic excitations are composed of a coherent su- 
perposition of many single-particle excitations. In our 
minimal model, we solve the eigenvalue equation (21), 
and the eigenvectors p™ (which depends on lu parametri- 
cally) describe how the transitions between noninteract- 
ing orbitals form the transitions in the interacting sys- 
tem. \p™\ is the percentage of a noninteracting tran- 
sition in the transition of the interacting system. Two 
typical cases are plotted in Figs. 6 and 7. 
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FIG. 6. Eigenvectors |p™| 2 of the first two excitonic tran- 
sitions for KP parameters a — 0.5, 6 = 0.5, Vo = 20. For 
BSE, A = 0.25 in Eq. (2). For TDDFT (soft-Coulomb), 
A 3C = 2.39 in Eq. (23). 
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FIG. 7. Eigenvector \p1 v \ 2 of a nonexcitonic excitation in the 
continuum part of the spectrum. The model parameters are 
the same as those in Fig. 6. 



Figure 6 clearly shows that excitons are collective ex- 
citations which are formed by mixing a wide distribution 
of single-particle transitions. As expected, the lowest ex- 
citon eigenfunction is nodcless and the second excitonic 
cigenfunction has a single node. With purely parabolic 
bands, the results from the Wannier model would be re- 
covered, as we will show below. In contrast, the tran- 
sitions in the continuum shown in Fig. 7 have a strong 
single-particle character (the two peaks arise from the ±fc 
degeneracy in the KP model). 

Equation (21) is equivalent to the SPA for finite sys- 
tems, which ignores the coupling between excitations and 
de-excitations (TDA). We also investigated what hap- 
pens when we do not make the TDA, i.e., when we work 
instead of Eq. (21) with the full equation for the two- 
band model, Eq. (19). As long as we describe relatively 
weakly bound excitons that are not too far below the 
band gap, we find that the difference between the two 
methods is very minor. 

However, we also discovered that, under rare circum- 
stances, Eq. (19) can lead to TDDFT excitonic bind- 
ing energies that are purely imaginary. In our mini- 
mal model, such instabilities arise when the interaction 
strength A in Eq. (23) and (24) increases so that the exci- 
tonic binding energy becomes greater than the band gap. 
This situation is comparable to the well-known triplet in- 
stability in TDDFT, for which the TDA generally leads 
to an overall better behavior; 31 for excitons binding ener- 
gies in our minimal model, we draw similar conclusions. 



B. Analysis of the coupling matrix 

The BSE scheme is commonly implemented within an 
adiabatic scheme (see Appendix A); as we have seen, it 
produces a series of bound excitons. Since we assumed 
that the KP model is the Kohn-Sham ground state in 
TDDFT and the GW quasiparticle ground state in BSE, 
the difference between TDDFT and BSE becomes easily 
comparable, since the central equation to be solved have 



the same form, Eq. (21). The if (*.#)(""*) coupling matrices 
for TDDFT and BSE are 



F 



(ij)(mn) _ 



TDDFT, k,k 



, = 2((ij\f H \mn)) + 2((ij\f xc \mn)), (25) 



ft 



BSE,k,k' 



2((u'|/h|™«)) - ((im\W\jn)}, 



(26) 



where /h is the Hartree kernel (the ID soft-Coulomb in- 
teraction), and W is the screened interaction. Aside from 
the change from / xc to W , the most prominent difference 
between BSE and TDDFT is the order of the indices for 
W in Eq. (26). Since the noninteracting ground-state 
wave functions have the Bloch form (11), we can see from 
Eq. (20) that F^^ mn ' has a strong k — k' dependence; 
in Fig. 8 this shows up as a dominance along the diago- 
nal. By contrast, this k — k' dependence is clearly absent 

in i^DDP"? ' as demonstrated in Fig. 9. 

((ij|/ xc |mn)) and ((im\W\jn)} with only vertical tran- 
sitions can be expressed in momentum space as 

mfxc\mn)) = - = °> G > G ') 



G,G' 



<j,k|e* G - r |*,k)(m,k' 



-iG'-r 



n,k'), (27) 



((im\W\jn)) = i ^ IF(q = k - k' + G , G, G') 

G,G' 

x (j,k e i(q+G) - r |n,k'\ /m,k' 



-*(q+G')-r 



i, k 



(28) 



The xc matrix (27) only depends on the long-range 
(q = 0) behavior of its momentum space representa- 
tion / xc (q, G, G'), while the W matrix (28) also depends 
on other q values in its momentum space representation 
W(q, G,G'). It is impossible to find an adiabatic f xc 
that reproduces the BSE coupling matrix as in Fig. 8, 
since W(q, G, G') has an extra degree of freedom over 
fxc(q = 0, G, G'). One can only hope to reproduce a por- 
tion of the BSE coupling matrix with adiabatic TDDFT 
(as pointed out in Rcf. 14), or make the xc kernel fre- 
quency dependent so that the information from the q- 
dependence in W(q, G, G') is mapped into the frequency 
dependence in f xc (q = 0, G, G',ui), 

Considering the nature of the objects involved in this 
mapping, a highly nontrivial frequency dependence in 
fxc is required to reproduce a series of bound excitons. 
For example, one can easily construct an / xc which re- 
produces a given series of bound excitons by using a dif- 
ferent contact kernel in the region — < uj < uii + rjf' 
surrounding each exciton at frequency uif. 

f xa (x,x',uj) = A l 8{x-x l )9[uj-(uj l ~ri~)}0{(uj l + ril)-uj]. 

(29) 

Here, the Aj's are parameters which are adjusted so that 
a TDDFT calculation with the frequency-independent 
kernel f xc (x,x') = Ai6(x — x') would produce Wj as the 
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-2 -1 





k 



with n° k k 



(o) 



,(0) 



-k -k ~ ' ano - F has the symmetry 
F k , k > = F_ k ,,_ k '= F k ,. k = F_ k , k , = F kt _ k ,. Within 
second order perturbation theory, there is at most one 
discrete eigenvalue of +F in the limit where k, k' be- 
come continuous. 42 Though this case does not correspond 
to the matrices that would occur in real calculations, it 
indicates the close relationship between the properties of 
the coupling matrix and excitonic effects. 

It is also possible to derive properties of the discrete 
eigenvalues if k and k' are completely decoupled in the xc 



(vc) (vc) 



F 



(vc) (vc) 



kernel. Owing to the symmetry F k k ; 
p(vc)(vc)* implied by Eq such 

separable kernels can 

only have the form 



F, 



k,k' 



±A(k)A*(k') 



(30) 



For an excitation below the band gap with frequency w, 
we can show 42 that it must satisfy 



FIG. 8. Contour plot of the coupling matrix j-fegg' 1 "^ | 
model parameters are the same as those in Fig. 6. 
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FIG. 9. Contour plot of the coupling matrix J^thdft 
model parameters are the same as those in Fig. 6. 



lowest excitonic binding energy. Such an / xc is of course 
completely ad hoc, but the fact that the excitonic se- 
ries can be reproduced in this way demonstrates that the 
inclusion of the frequency dependence would greatly im- 
prove the flexibility of the TDDFT scheme. 

On the other hand, within the adiabatic approxima- 
tion the characteristics of the i^xc coupling matrix 
are important for excitonic effects. To emphasize this, 
we now show that in very special cases the number of 
discrete excitonic eigenvalues can be derived. Consider a 
real matrix + F, where f2° is a diagonal real matrix 



k * 



(31) 



where the sum is carried out over the first Brillouin zone 
(FBZ). Equation (31) shows that Eq. (30) must have the 
negative sign in order to have bound excitons. The left- 
hand side of Eq. (31) is monotonically increasing with oj, 
so for separable kernels of the form of Eq. (30), there is 
only one bound excitonic solution. 

As shown in Fig. 9, TDDFT coupling matrices lack the 
strong dependence oik — k' as in BSE coupling matrices. 
Expanding the TDDFT coupling matrices into a power 
series of separable matrices and truncating at the first 
order would be a reasonable approximation, explaining 
why TDDFT produces fewer bound excitons (if any at 
all) than many-body methods such as BSE. 



C. Dimensionality considerations 

The contact xc kernel and the soft-Coulomb xc kernel 
in Eqs. (23) and (24) have the following simple form in 
momentum space: 



/|°(g, G, G') = -2A SC K (y^\q + G|) 6 G ,G>, 
f™\q,G,G') = -A 



(32) 



where q G FBZ, G and G' are reciprocal lattice vectors, 
and K is the modified Bcsscl function of the second kind. 
It is customary to refer to the matrix elements where 
G = G' = as 'head', G = or G' — as 'wings', and 
G^0,G'^0as 'body'. 

The 3D Coulomb potential has the form 4ir/q 2 in mo- 
mentum space. However, in ID systems there is no real 
Coulomb interaction which behaves as q~ 2 for q — > 0, 
and one has to use the soft-Coulomb interaction instead. 
Though there are many flavors of the soft-Coulomb inter- 
action, they all have the same logq behavior for q — > 0. 
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However, the linear response function x always behaves 
as q~ 2 for q — > and does not depend on the dimen- 
sionality. This renders quantities like the macroscopic 
dielectric function (12) ill-defined for strictly ID systems. 
Furthermore, the bootstrap xc kernel 20 and other xc ker- 
nels that depend on the cancellation of the 3D Coulomb 
q~ 2 singularity will not work as designed in strictly ID 
and 2D systems. Therefore Im(ejvf) shown in Fig. 4 and 
5 are calculated at a small but finite q. 

The coupling matrix Fhxc l ' c ' can be written in mo- 
mentum space as 



R 



(vc) (vc) 



HXC.fci' 



\ E M« = °) 5 G,G' + fxc(q = 0, G, G')] 



fi 



G,G' 



(c,k\e lGx \v,k) (v,k' e~ lG ' x c,k'\. (33) 



For any xc kernel that behaves as q~ 2 for q 0, one 
can further simplify the calculation by ignoring the so- 
called local field effects, 43 i.e. instead of summing over 
G and G' in Eq. (33), only the head is considered. In 3D 
systems, a prominent example is the long-range kernel 
—a/q 2 , which is obtained as an effective xc kernel with 
only head matrix elements from inverting the BSE of 
contact excitons. 15 

On the other hand, any xc kernel that diverges more 
slowly than q~ 2 for q — > changes the spectrum only 
through the local field effects, i.e. all G and G' must 
be summed in Eq. (33). In other words, effective xc 
kernels with only the head are not feasible in strictly 
ID systems due to the asymptotic behavior of the soft- 
Coulomb potential discussed above. For ID systems with 
G = 0. we have 



i. ki 



<3->0 



Oiq 1 ), 



v G =o(q) V ~° O(logg), 



(34) 



and /xc's with ID long-range behavior such as the soft- 
Coulomb kernel also behave as O(logg). Considering 
Eq. (33), these asymptotic properties imply that the 
head and wing contributions to _FW)( m ") always vanish 
in strictly ID systems for physically meaningful xc ker- 
nels. Due to these dimensionality restrictions, the xc 
kernel changes the strictly ID system only through the 
local field effects. 

In 3D the head contribution to the coupling matrix 
F H xc is much more important than the local field effects, 
which is the reason that long-range kernels (with nonzero 
head) work much better than local xc kernels (with van- 
ishing head) such as ALDA. In our strictly ID model 
system, the head contribution is zero even for the BSE, 
and thus the long-range kernel does not outperform local 
kernels such as the contact kernel. 

These peculiarities only occur when one considers 
strictly ID and 2D systems. In a more realistic pic- 
ture, one encounters quasi-2D systems 44 and quasi-lD 
systems (such as quantum wires with finite radius or 



nanotubes 45 ), in which the movement of electrons is con- 
fined in certain directions such that the transverse mo- 
tion can be averaged in comparison with the longitudi- 
nal motion. Though these systems show low-dimensional 
characteristics in various properties due to confinement, 
in the limit of q — > they eventually differ from strictly 
low-dimensional systems. 



IV. THE WANNIER MODEL IN TDDFT 

Our minimal model and the Wannier exciton picture 
can be connected by considering the Fourier transform of 
Eq. (21). We define an effective two-body potential V e h 



via the Fourier transform of _F! 



(vc)(vc) 



y oh (R,R') = 



2tt 



Hxc.k.k' 



-ik-Rp("c)(oc) „ik'-R' 



Hxc,k,k 



,e 



k,k'£FBZ 



2m r 



(35) 

where R is a direct lattice vector. The Fourier transform 
of the density matrix is 

p(R,u) = p cv (B.,u) ^e- !k 'V( W ). (36) 



Since Wannier excitons extends over many lattice con- 
stants, we approximate R as a continuous variable r. 
Assuming the effective mass approximation, Eq. (21) 
becomes 

p(r,u) + / dV V r eh(r J r / ,w)p(r , ,w) = Ep(r,u), 

(37) 

where E is the excitonic binding energy, and the integra- 
tion is carried out over all space. We call Eq. (37) the 
TDDFT Wannier equation, since it has the same form 
as Eq. (1). With proper choice of the approximated 
xc kernel, the nonlocal effective electron-hole interaction 
potential 14h supports bound excitonic states. 

Since the BSE and the TDDFT are formally similar 
within the minimal model, Eq. (21) can also be applied 
to the BSE results. Fig. 10 shows the effective interaction 
potential V ch for TDDFT and BSE. 

The TDDFT Wannier equation provides an intuitive 
way of describing the effective nonlocal electron-hole in- 
teraction, and of explaining why adiabatic TDDFT usu- 
ally has fewer excitons than BSE and the Wannier model. 
However, in most cases the TDDFT Wannier equation is 
not suitable for quantitative use due to the approxima- 
tions involved. The approximation where we take the 
lattice vector R as a continuous variable assumes that 
the exciton radius is much larger than the lattice con- 
stant; this works fine in most cases we tested. But the 
effective mass approximation where w q is approximated 
by q 2 /2m r is only good for transitions near the band gap, 
thus requiring these transitions of the noninteracting sys- 
tem to dominate the exciton, which is equivalent to the 
exciton extending over many lattice constants. One ob- 
tains the — V 2 /2m r term in Eq. (37) from q 2 /2m r in 



10 



(a) 



V o 




(b) 



I — — 1 u> 

6 d q a 
■ODD 

4 3 a a 
■ODD 

2 3 an 



-6 



- 1 — cp — 1 — ax - ' — 

a a o 4 

D 
a a a e 

D 
n) a a„e 



O D D/D7j i; 



3 , 
4 ] 

ODD 

3 o a 

m, m, - 



G a o * 

» D O D. 

a a o * 



0.4 
0.2 



-0.2 
-0.4 
-0.6 
-0.8 
-1 




x 



FIG. 10. Contour plots of V ch (r,r') for (a) BSE, and (b) 
TDDFT with the soft-Coulomb kernel. The KP system is the 
same as in Fig. 6. 



the limit where the lattice constant a + b — > 0, and this 
approximation is not valid for most systems. 

Although Vch is a nonlocal potential, in most cases 
we find that the V^h's for both TDDFT and BSE are 
dominated by the diagonal part, so the exciton problem 
is in analogy to one-body systems. Fig. 11 shows the 
diagonal part of V e h, which can be taken as the effec- 
tive one-body potential. The Wannier model in ID has 
the soft-Coulomb potential, which supports an infinite 
number of bound excitons (the soft-Coulomb interaction 
is fitted so that the binding energy of the first exciton 
matches that of the BSE). We find in general that the di- 
agonal parts of Kh for both BSE and TDDFT are much 
more shallow than the soft-Coulomb potential, and the 
TDDFT one is more narrow than the BSE one. Thus, 
BSE and TDDFT are not able to produce a complete ex- 
citonic Rydberg series, and TDDFT in general produces 
fewer bound excitons than the BSE. 

The TDDFT Wannier equation is not suitable for 
quantitative use for most of our model systems, de- 
spite the success of the Wannier model in describing real 
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FIG. 11. Diagonal part V a h(r, r) of those shown in Fig. 10. 



semiconductors. 46 ' 47 Since the approximations involved 
in Eq. (37) require that the exciton radius is large com- 
pared to the lattice constant, this suggests that this dis- 
crepancy is due to the special nature of ID systems: 
namely, for similar effective masses the exciton radius 
in ID is much smaller than in 3D and 2D. 27 



V. CONCLUSION 

The purpose of this paper was to construct a trans- 
parent and accessible minimalist model system that pro- 
duces excitonic effects in a non- ad hoc fashion using 
TDDFT. The model, as presented here, is not intended 
to be a testing ground for xc kernels. Thus despite the 
dimensionality restrictions for the strictly ID system, our 
results in Sec. Ill carry over to 3D systems. 

With our minimal model, we show that adiabatic 
TDDFT is capable of producing bound excitons through 
the local field effect even when the xc kernel is local 
in space, provided the strength of the kernel is strong 
enough. This statement is still true in 3D; however, due 
to the non- vanishing head contribution of the exact / xc , 
we expect that the deviation of the effective interaction 
strength of a local / xc from the real, nonlocal / xc be- 
comes larger than the in our strictly ID model. In this 
sense the long-range kernel, though very favorable, is not 
a necessary condition for excitonic effects. 

We show the connection between TDDFT and the 
Wannier model for excitons by deriving the TDDFT 
Wannier equation, which describes a real-space system 
featuring a nonlocal effective electron-hole interaction. 
Such a connection intuitively demonstrates how adiabatic 
TDDFT generally produces fewer bound excitons than 
BSE, and does not have a complete Rydberg series. The 
eigenvectors of the excitonic excitations in the minimal 
model clearly show their collective nature, which is not 
obvious from the Wannier model alone. Excitonic in- 
stabilities may show up in TDDFT with approximate xc 
kernels, and this suggests that the TDA tends to be more 
reliable for excitons than the formally exact method. 



11 



The frequency dependence of the exact xc kernel, 
/ xc (r, t',u), is usually ignored. Despite the fact that 
adiabatic xc kernels have met with considerable recent 
success in producing optical spectra of insulators and 
semiconductors (see the discussion in the Introduction), 
they are incapable of producing excitonic Rydberg series. 
Our model system gives an explanation for why this is 
the case. This failure of the adiabatic approximation 
for /xc is quite different from that which is responsible 
for the inability of adiabatic TDDFT to produce double 
excitations in finite systems or certain classes of charge- 
transfer excitations. 48 ' 49 This calls for continuing efforts 
in the search for nonadiabatic xc kernels for excitons. 



ACKNOWLEDGEMENT 

This work was supported by NSF Grant No. DMR- 
1005651. 



Appendix A: The Bethe-Salpeter equation 

Electrons and holes near the Fermi surface are well 
described in the quasiparticle picture. The quasiparticle 
Green's function G is related to that of the noninteracting 
system, Go, through the use of the self-energy S: 

G(12) = G (12) + J d(34) Go(13)E(34)G(42), (Al) 

where the arguments denote sets of space and time vari- 
ables. A widely used approximation for the self-energy is 
the GW approximation: 1,2 

E(12) = iG(12)W(12), (A2) 

where W is the screened interaction, 

W{r,r',oj)= I d 3 r" e^fcr" »v(r" >'), ( A 3) 



v is the bare Coulomb interaction, and the inverse dielec- 
tric function e _1 is obtained as 

e- 1 {v,v',iu) = d{r-r')+ f dV v(r, r")x(r", r', u). 



(A4) 

The linear response function \ can be calculated by its 
Lchmann representation: 

, y,;(r)^(r)^(rO^(Q 
X (r, r , u = 2^ tt; FT^~ wi ~ & ' A5 

where ipi are quasiparticle states, Ei arc quasiparticle en- 
ergies, and ft are occupation numbers. In practice, eval- 
uating \ through Eq. (A5) can be quite time-consuming, 
and x is often calculated with the plasmon-pole model. 50 
However, our minimal model is simple enough to allow 
us to use Eq. (A5) directly. 



The GW quasiparticle Green's function obtained from 
Eq. (A2) misses important dynamical many-body effects, 
such as the electron-hole (excitonic) interaction. The 
two-particle Green's function includes these effects. The 
BSE 3 ' 4 describes the relation between the four-point po- 
larization function £(1234) of an interacting system and 
the corresponding object of the quasiparticle system: 

£(1234) = £ (1234) + J d(5678) £ (1256) 

xA(5678)£(7834), (A6) 

in which Lq is 

£ (1234) = iG(13)G(42) (A7) 

and assuming the GW approximation, the kernel K is 

#(1234) = <5(12)<5(34)v(13) - <5(13)<5(24)W(12). (A8) 

Here, v denotes the Coulomb interaction with the long- 
range part removed. 4 In practice the BSE is often solved 
in the transition space, which is spanned by single- 
particle excitations. A four-point function such as £ then 
becomes 

£ fo)(lml) (w) = j dan . . . x 4 £(r 1 r 2 r 3 r 4 ; w) 

x ^(r^^WMM^), (A9) 

where the <p's can be any complete basis set. Eq. (A6) 
in the transition space becomes 

L(«X™0( W ) = [h cxc (uj) Iw]^ (mn) (/ m - /»), (A10) 

where the excitonic Hamiltonian matrix is 

H^ mn \u) = (Ej -Ei- to)6 im S jn 

+ (fi- fj)K^ mn Huj). (All) 

We make the adiabatic approximation for H^ic and 
arrive at the following eigenvalue problem: 

van 

and £w)( m ") can be expressed in terms of these eigen- 
vectors by 



w^ xc — 



(A13) 



Only the transitions between the valence and conduction 
bands contribute. Within our two-band model, the exci- 
tonic Hamiltonian has the following block matrix form: 



i£ 



E c - E v + K ( - VC X VC '> \ 

_j((vc)(cv)* ^ — E c — X( vc )( vc )* J ' 

(A14) 
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Ignoring the off-diagonal part in Eq. (A14) is equivalent 
to the TDA. 

As shown in Sec. Ill, it is possible that instabilities 
show up in the full BSE results when the underlying 
ground-state calculation is not exact. Such instabilities 
in the minimal model are an artifact originating from the 
assumption that the solution of the KP model constitutes 
the ground state of the many-body system. However, this 
is not a matter of great concern in practice. 

In principle, the transition space spans all possible 
combinations of valence and conduction orbitals, includ- 
ing nonvertical transitions connecting different Bloch 
wavevectors. The kernel K = v — W of the BSE in mo- 
mentum space, Eq. (A8), has the following ingredients: 



v (v)(mn } = 1 £ UG (q)£ q>k ._ ki+Go * q)k „_ km+Go 



n 

i(q+G)-r| ; k .\/ mki 



-»(q+G)-r 



TL 7 k n / 



(A15) 



and 



W immn) = _J_ £ W G , G , (q)«J q , kj - k „ +Go <5 q , ki - km+Go 



G,G' 



x (j,k,| e ^ +G )- r |n,k„)(m,k r , 



,-i(q+G')T 



i, ki 
(A16) 



Only the excitations with the same momentum transfer 
q are coupled due to the S functions in Eq. (A15) and 
(A16), so we only need to include vertical transitions in 
the calculations for optical properties. 
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